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1. Introduction 

Lattice QCD simulations are routinely done 
nowadays using highly- improved actions, which 
are designed to remove the leading errors arising 
from the discretization of the continuum theory. 
Improved actions can be designed from both per- 
turbative and nonperturbative considerations. In 
this review I describe techniques for doing the 
higher-order perturbation theory (PT) calcula- 
tions that are necessary in the design of highly- 
improved lattice discretizations for gluons, stag- 
gered quarks, and heavy quarks. Recent results 
obtained with these methods are also reviewed. 

Much of the work described here is part of the 
program of the HPQCD collaboration, a major 
goal of which is to make precision calculations of 
hadronic matrix elements relevant to 6-physics ex- 
periments. In order to fully realize the potential 
impact of these experiments on the parameters of 
flavor mixing requires calculations of the relevant 
hadronic matrix elements to a few percent accu- 
racy. The CLEO-c program also presents an enor- 
mous opportunity to validate lattice QCD meth- 
ods for b physics, by testing predictions for anal- 
ogous quantities in the charm system. 

These stringent requirements for accuracy and 
timeliness are unlikely to be met without signif- 
icant algorithmic developments, in both the ef- 
ficiency of unquenched simulations, and in the 
technical challenges posed by lattice PT. The 
development of an improved action for stag- 
gered quarks |1I2| has at last made accurate un- 
quenched simulations feasible |3] at dynamical 
quark masses that are small enough to allow for 
reliable chiral extrapolations to the physical re- 
gion 4 . This is one particularly striking success 
of the perturbative analysis of lattice discretiza- 



tions. More generally one must do perturbative 
matching calculations for a wide array of coupling 
constants, action parameters and hadronic matrix 
elements. 

The scope of the charge to lattice PT is set 
by two expansion parameters: oAqcd, where a 
is the lattice spacing and Aqcd is a typical low- 
energy scale; and the strong coupling a s (q*), eval- 
uated near the ultraviolet cutoff (q* ~ 1/a) at 
which the lattice theories are to be matched onto 
continuum QCD. Affordable unquenched simula- 
tions can only be done for lattice spacings around 
0.1 fm, where these two expansion parameters 
have about the same value: 

a s (1/a) » aAqco w 0.2-0.3 [a w 0.1 fm]. (1) 

Hence to reduce systematic errors to a few percent 
requires lattice discretizations that are accurate 
through 0(a 2 ), which is trivial to do by numerical 
analysis, and matching of the resulting action and 
operators must be done through 0(a 2 ). 

The latter requirement is highly nontrivial, and 
in fact only a few two-loop PT calculations in lat- 
tice QCD have ever been done. In many ways per- 
turbative calculations with a lattice regulator are 
much more difficult than with dimensional regu- 
larization: the Feynman rules for lattice theories 
are exceedingly complicated, even for the simplest 
discretizations, with a tower of contact interac- 
tions that leads to a proliferation of Feynman di- 
agrams that are not present in the continuum. 

One has the challenge of doing many two-loop 
matching calculations, and doing these for the 
most complicated actions yet designed. Moreover 
one may anticipate further evolution in the ac- 
tions that are actually used in simulations, as in- 
vestigations continue into the optimal discretiza- 
tions, hence one will want to be able to redo many 
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perturbative calculations as the actions evolve. 

Fortunately there exists several techniques 
which make higher-order lattice PT more man- 
ageable. In this review I survey two very differ- 
ent and somewhat complementary methods: one 
for conventional Feynman diagram evaluation of 
loop corrections, and another approach in which 
one does Monte Carlo simulations of the lattice 
path integral in the weak-coupling regime. 

The linchpin of our program for diagrammatic 
PT is the automatic generation of the Feynman 
rules for lattice actions. This method was de- 
veloped long ago by Liischer and Weisz 0, and 
is described in Sect. 2. What is new in our 
work is that we have aggressively applied this 
method to higher-order calculations for highly- 
improved actions for gluons, staggered quarks, 
and heavy quarks. This is in order to analyze 
the MILC simulations of (2+l)-dynamical quark 
flavors [3], which were done for the 0(a s a 2 )- 
accurate Symanzik gluon action, and the tree- 
level 0(a 2 )-accurate staggered quark action due 
to Lepage pQ, both with tadpole improvement. 

A major effort of the past year has been de- 
voted to a complete third-order determination of 
a s (Mz), the first from lattice QCD, which should 
provide one of the most accurate determinations 
of the strong coupling. This requires a set of two- 
and three-loop calculations, which are described 
in Sects. 3 and 4. 

Some other recent PT calculations are briefly 
reviewed in Sect. 5; a comprehensive review of re- 
cent work cannot be done in this short space, so I 
have mainly restricted this summary to work done 
with automated methods. Monte Carlo meth- 
ods for the extraction of perturbative expansions 
are described in Sect. 6. Some conclusions and 
prospects for the future are found in Sect. 7. 

2. Automatic vertex generation 

Lattice Feynman rules are vastly more compli- 
cated than in the continuum, even for the sim- 
plest discretizations (the continuum four-gluon 
vertex for instance has six terms, while for the 
Wilson gluon action the expression spans dozens 
of lines). The complexity of these rules grows ex- 
tremely rapidly even for modest increases in the 



complexity of the action, and with the number of 
lines in the vertex. A path in the action with £ 
links will generate a vertex function for r gluon 
lines with a number of terms bounded by |S] 

n r>i = 21(1 + l)...(i + r- l)/(r - 1)!. (2) 

The C(a 2 )-accurate gluon action has I — 6, while 
for third-order expansions one requires vertex 
functions with r — 6 gluons, which have up to 
"6,6 = 5544 terms (cancellations lead in practice 
to expressions of about half that size). 

On the other hand remarkably simple algo- 
rithms can be developed to automate the gen- 
eration of the Feynman rules, for essentially arbi- 
trarily complicated lattice actions. The method 
which we use combines algorithms of Liischer and 
Weisz and Morningstar (i . Other algorithms 
have been developed by Panagopoulos and col- 
laborators 0, and Capitani and Rossi 

The algorithm has a very user-friendly inter- 
face. One need only specify the action in an ab- 
stract form, according to the link paths that it 
contains. The algorithm then Taylor expands the 
link variables, collects terms of the desired order 
in the coupling, keeping track of Lorentz and color 
indices, and Fourier-transforms the fields. A final 
expression for the vertex function is printed in a 
language that is suitable for use in routines where 
the Feynman diagrams are coded, from a combi- 
nation of vertex functions, allowing for instance 
a numerical integration over the loop momenta. 
(Automatic generation of the Feynman diagrams 
themselves can also be done 9 .) 

To turn this prose into a computer program 
that can handle a very complicated action, it is 
helpful to observe that one can build up the Feyn- 
man rules from a convolution of the rules for sim- 
pler elements .6.. To this end, let's consider the 
world's simplest "action," consisting of a link in 
a single direction, summed over lattice sites 

X X 

We can immediately write down the Feynman 
rules for this "action," since we need only to take 
the Taylor series expansion of the lone exponen- 
tial. The vertex function V for r gluons is 
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Kmsymm ({&ti Mi: a i}) ~ ^ X 

(2tt) 4 <5 (fci + fc 2 + • • • + K = k tot ) x 

4<KAi = M2 = • • ■ = Mr = i) x {T ai T a ' 2 . . . T a "} 
r! 

x exp[i(i/ci • at + \k 2 ■ at + . . . + \k T ■ at)], (4) 

where we first compute an unsymmetrized vertex, 
with the gluon labels assigned in a fixed order to 
the link matrices (the color trace has not yet been 
taken); a sum over permutations of the labels 
is done once the unsymmetrized amplitude for 
the complete action has been constructed. The 
origin of each term in Eq. Q is simple: mo- 
mentum conservation (allowing for momentum 
transfer k to t into the vertex), a 1/H coming from 
the Taylor expansion, and phase factors from the 
Fourier transforms of the gauge fields. The spin 
(5-function arises because all gluons must be po- 
larized along the direction of the single link. 

To build the vertex function for a more complex 
action one can convolute the vertices of individ- 
ual links. This is illustrated by another simple 
example, in Fig. ^ which shows the generation of 
the unsymmetrized vertex function for two gluons 
(r = 2) for an "action" with two links (£ = 2). 
One simply assigns an ordered set of gluon labels 
to the ordered set of links (so as to respect their 
non-commutativity), in all possible ways. 

Gluon actions with relatively short link paths 
do not require a program for convoluting the rules 
of individual links; this can easily be done analyt- 
ically, and Liischer and Weisz provide an explicit 
realization of an optimal algorithm in that case. 
[They also discuss how to automatically gener- 
ate Feynman rules for gauge-fixing, although for 
practical purposes this can be done by hand.] In 
more complex cases however it is advantageous 
to explicitly code the convolution theorem, ap- 
plied to basic operators in the action. Consider 
for instance the NRQCD action 

Wd =*»(l -!«,)(! + _) U\ 

1 + ^Mj (l-H^V, (5) 
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Figure 1. Generation of an unsymmetrized vertex 
function by convolution. 



where SH has many operators 

a6H=—^a--(AxE-Ex A)— ^tr-B+. . . .(6) 

8M 2 y ' 2M y ' 

The vertex functions for this extremely complex 
action can easily be built-up by convolution of the 
rules for simple operators, such as A^ 2 ) and E and 
B 6 . We have done a number of calculations for 
NRQCD using exactly this procedure. 

3. Two-loop renormalized coupling 

The convergence of perturbative expansions is 
greatly improved by using a renormalized cou- 
pling 10 , such as defined by the static-quark po- 
tential, according to 

V(q) = ~4TrC F a v (q)/q 2 . (7) 

We have computed the connection between ay (q) 
and the bare lattice coupling ao for improved ac- 
tions in two steps We first follow Liischer 
and Weisz ^21 > wn0 used a background-field tech- 
nique to compute the relation between ao and 
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a MS f° r ^ ne Wilson gluon action (the matching 
for Wilson and clover fermions was done in 0). 
We then use the relation between a^g and ay, 
which is known through third-order |13) . 

As in the continuum, background-field quanti- 
zation reduces the number of independent renor- 
malization constants. One computes background- 
and quantum-field two-point functions on the lat- 
tice, r B and T A respectively, where 

T B (q,-q)^ = -S ab 3f [l~v{q)]/gl 

r A (q, -q)% = -5 ab q 2 [3 (1 - u(q)) + A„] , (8) 

where go and Ao are the bare lattice coupling and 
gauge parameter. Analogous renormalized quan- 
tities <7f^g, vr, etc. are defined in the MS scheme. 
The couplings in the two schemes are related by 

fl|g = [l-i'fl(5)]/[l-K?)]x 9 g. (9) 

To solve this equation one must account for the 
implicit dependence on the couplings induced by 
the relation \r = [1 — w^(g)]/[f — u{q)\ x A . 

At two-loop order one must evaluate 3f pure- 
gauge diagrams, and 18 diagrams with internal 
fermion lines. In addition one must compute a 
number of one-loop diagrams that are induced by 
the 0{g 2 ) expansion of the tadpole renormaliza- 
tion factors in the gluon and quark actions. 

We evaluate loop integrals by numerical inte- 
gration using VEGAS. This contrasts with |12I7| . 
where integrals for unimproved actions were re- 
duced analytically to a small set of primitive inte- 
grals. In the present case the much more complex 
vertices make analytical treatments problematic 
(see however |14p: moreover there is no gauge in 
which the improved gluon propagator is diago- 
nal in its Lorentz indices, which further compli- 
cates analytical integration. We find that numer- 
ical integration gives results of sufficient quality 
in reasonable computer time. We do integrations 
at several values of aq, and extrapolate to the 
continuum limit, as illustrated in Fig. [21 We re- 
quire the fit errors in third-order coefficients to be 
smaller than the systematic error from the uncal- 
culated fourth-order corrections, for the couplings 
relevant to simulations. 

We also do calculations for several choices of 
the gauge parameter, to explicitly verify gauge- 
invariance of the final results. This requires the 
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Fit = -2.6143(11) 

L&W = -2.6126215... 
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Figure 2. Two-loop results for v{q) (Wilson glue). 



MS matching function vr at two loops for ar- 
bitrary gauge parameter. This was done for the 
pure-gauge theory by Ellis |15) , and we have done 
the two-loop fcrmionic part for arbitrary A_r (this 
was previously known only in Feynman gauge jxj ) • 
The relation between the couplings is given by 

"o = a v (q) [l - v 1 (q)a v (q) - v 2 (q)a v (q)] , (10) 

plus corrections of 0(a.y), where 

Vx(q) = (A)/4tt) ln(ir/aq) 2 +«i,o, 

v 2 (q) = (/3i/16^ 2 )ln(7r/ag) 2 - v\{q 2 ) + v 2 . . (11) 

For the improved actions used by MILC we find 

wi,o = 3.57123(17) - 1.196(53) x Kr 4 ^/, 

u 2 ,o = 5.382(39) - 1.0511(51)^/. (12) 

4. Wilson loops to third order 

One can extract the value of the strong cou- 
pling from lattice simulations of short-distance 
observables. To this end we have evaluated the 
three- loop PT expansion of Wilson loops ex- 
tending a clever approach in in order to reduce 
the number of Feynman diagrams. We evaluate 
a vacuum-to-vacuum transition amplitude, in the 
presence of an external "current" generated by 
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the Wilson loop itself 

(W R , T ) 



- P W R: 



[atp = 0], (13) 



using our vertex generators to make the rules for 
the "extended" action Siat+pW^T • This reduces 
the number of three-loop diagrams by more than 
half compared to a "direct" evaluation of (Wr,t)- 
The third-order expansion has 15 gluonic three- 
loop diagrams, and 19 fermionic ones. In addition 
there are a number of one- and two-loop mean- 
field counterterm diagrams. We quote results for 

— lnW R:T = Woa v (q*)[l + ria v +r 2 cty}+. . . .(14) 

The gluonic parts r^ g are given in Tabled which 
demonstrates convergence of the renormalized PT 
expansion through three loops. 
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Figure 3. Quenched coupling vs. lattice spacing 



Table 1 

Perturbative expansions for small Wilson loops 
for the 0(a s a 2 )-improved gluon action. 



i?xT 


aq* 


w 






lxl 


3.33 


3.0684 


-0.7753(2) 


-0.722(39) 


1x2 


3.00 


5.5512 


-0.6202(4) 


-0.407(40) 


1x3 


2.93 


7.8765 


-0.5335(8) 


-0.245(44) 


2x2 


2.58 


9.1998 


-0.4934(10) 


-0.030(51) 



Results for the quenched coupling extracted 
from simulations at different orders in the pertur- 
bative expansion are shown in Fig. [3| 16 . This 
indicates convergence through fourth order, with 
the results from two gluon actions at nth-order 
differing by 0(1) x a^ +1 , all the way from n = 1 
through n = 3. We are currently analyzing the 
results for the MILC (2+l)-flavor configurations. 

5. Other recent results 

The 0(a 2 )-improved staggered-quark action re- 
moves taste-changing interactions at tree level 
by suppressing the coupling of quarks to high- 
momentum gluons |1I2| . The ancient problem of 
a bad perturbation theory for staggered quarks 
has also been shown to be thereby eliminated |17| . 
These effects can be removed at higher orders by 
introducing four-fermion counterterms, cf. Fig. 



Hascnfratz has found that taste-symmetry 
breaking can be further reduced by additional 
smearing and reunitarization of the gauge links in 
the staggered-quark action JHj. We have found 
that this improvement is also largely perturbative 
in origin |19| . Taste-symmetry violations were 
measured from pion splittings in quenched sim- 
ulations of staggered actions with different link 
prescriptions, and calculations of the one-loop 
four-fermion counterterms were also done in each 
case. Perturbation theory correctly accounts for 
changes in the pion splittings with changes in the 
action, and the same optimal action is found per- 
turbatively and nonperturbatively. Algorithms 
for doing unquenched simulations with these ac- 
tions are also under investigation |20) . 

Another successful PT description of nonper- 
turbative simulation data comes from a computa- 
tion of the one-loop renormalized anisotropy for 
improved gluon actions [21]. Perturbative and 
nonperturbative determinations of the anisotropy 
were compared over a wide range of lattice spac- 
ings and anisotropics, the differences in all cases 
consistent with an 0(1) x oiy correction. 

An important application of lattice PT is 
matching of the Fermilab and NRQCD heavy- 
quark actions and associated currents [221231241 
23 . The first PT matching for the heavy-quark 
clover interaction was done in |26|. and new re- 
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Figure 4. Taste-changing effects for improved- 
staggered quarks at tree-level and one-loop. 



shown in Fig. 03 the intercept of the graph is the 
leading order coefficient co, while the slope of the 
curve resolves ci, and its curvature resolves C2. 
We subsequently did the three-loop expansions 
by Feynman diagrams , and the results are in 
excellent agreement, see Table EJ 
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suits were presented here by Kayaba . This is 
the subject of a concerted effort by the HPQCD 
collaboration; complete one-loop results for the 
action parameters are coming soon |23| . We are 
also well on the road to the two-loop kinetic mass 
for light and heavy quarks ^1] (see also |28|). 

6. Monte Carlo methods 

There is an attractive alternative to Feynman 
diagram analysis, which has largely been under- 
exploited. This is to "do" PT by doing conven- 
tional Monte Carlo evaluations of the lattice path 
integral, but in the "unconventional" weak cou- 
pling phase of the theory |29| . 

One simulates an observable, 

(0) = J[dU ll {x)][d^{x)dil;{x)]Oe'l 3S ^\ (15) 

over a range of weak couplings, say (3 w 9-60, in 
a finite (Planck!) box, where all lattice momenta 
are very large (except for possible zero modes, 
which can be eliminated by an appropriate choice 
of boundary conditions). One then fits the results 
to the series (O) = J^n c « a y (?*)• Third-order 
expansions of Wilson loops and the static-quark 
self-energy were obtained in quenched simulations 
of the plaquette action in |3U|. Some results are 



Figure 5. Monte Carlo data for a 5x5 Wilson loop, 
with Ki = -ln(VF 5x5 )/(20a(g*)) [Wilson glue]. 



Table 2 

Monte Carlo (MC) and PT results for r 2 , g (Wil- 
son gluon action). 



Loop 


MC 


PT 


1 x 2 


-1.34(8) 


-1.31(1) 


1 x 3 


-1.17(9) 


-1.13(5) 


1 x 4 


-1.04(9) 


-1.00(8) 


1 x 5 


-0.98(11) 


-1.05(12) 


2x2 


-0.71(8) 


-0.71(5) 


2x3 


-0.44(9) 


-0.48(9) 


3x3 


-0.12(9) 


0.15(17) 



An alternative method has been developed by 
the Parma group, using an explicit expansion of 
the Langevin equations in the bare coupling; they 
have recently obtained the unquenched third- 
order static-quark self-energy [2]. We have 
also recently extracted PT expansions from un- 
quenched weak-coupling simulations |32|. 
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7. Summary and outlook 

Automatic lattice perturbation theory meth- 
ods are remarkably simple and powerful. These 
have been used to do a number of higher-order 
calculations for highly- improved actions. We 
have recently done the PT for a third-order de- 
termination of a^a(Mz) from the MILC simu- 
lations of (2+l)-flavors of dynamical staggered 
quarks. Conventional Monte Carlo simulations 
in the weak-coupling regime may offer a simple 
alternative to diagrammatic PT, which has so far 
been under-exploited. Our primary goal is a two- 
loop determination of heavy-flavor physics. This 
is an ambitious program but the technology has 
been proven at the requisite order, and one can 
reasonably expect that much of this work will be 
done in the next few years. 

I am indebted to Peter Lepage, who has pro- 
vided much of the rationale for this physics pro- 
gram. I have done most of my PT calculations 
with Qucntin Mason, Matthew Nobes, and Peter 
Lepage. I have benefited from many discussions 
with Christine Davies, Junko Shigemitsu, Paul 
Mackenzie, Andreas Kronfeld, Aida El-Khadra, 
Richard Woloshyn, and Ron Horgan. 
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